%% setup
N_corr_sim = 100;% # simulation draws for the common shock
profit_cutoff_corr = profit_cutoff;

%% sample pairs of products
all_prod_pair_ind = nchoosek(1:N_dev, 2);
N_prod_pair_ind = size(all_prod_pair_ind, 1);
N_dev2 = min(N_dev2, N_prod_pair_ind);
ind_est_corr = randsample(N_prod_pair_ind, N_dev2);
dev2_ind = all_prod_pair_ind(ind_est_corr, :);

%% find their profits
dev_profit_all_dev2_f1 = reshape(dev_profit_all(:, dev2_ind(:, 1))', [N_mkt*N_dev2 1]);
dev_profit_all_dev2_f2 = reshape(dev_profit_all(:, dev2_ind(:, 2))', [N_mkt*N_dev2 1]);
dev_profit_one_dev2_f1 = reshape(dev_profit_one(:, dev2_ind(:, 1))', [N_mkt*N_dev2 1]);
dev_profit_one_dev2_f2 = reshape(dev_profit_one(:, dev2_ind(:, 2))', [N_mkt*N_dev2 1]);

dev_profit_all_dev2_f1_boot = reshape(permute(...
    dev_profit_all_boot0(:, dev2_ind(:, 1), :), [2 1 3]), ...
    [N_mkt*N_dev2 N_boot]);
dev_profit_all_dev2_f2_boot = reshape(permute(...
    dev_profit_all_boot0(:, dev2_ind(:, 2), :), [2 1 3]), ...
    [N_mkt*N_dev2 N_boot]);
dev_profit_one_dev2_f1_boot = reshape(permute(...
    dev_profit_one_boot0(:, dev2_ind(:, 1), :), [2 1 3]), ...
    [N_mkt*N_dev2 N_boot]);
dev_profit_one_dev2_f2_boot = reshape(permute(...
    dev_profit_one_boot0(:, dev2_ind(:, 2), :), [2 1 3]), ...
    [N_mkt*N_dev2 N_boot]);

%% construct insturments
fc_iv0_dev2_f1 = reshape(permute(fc_iv0(:, dev2_ind(:, 1), :), [2 1 3]),...
    [N_dev2*N_mkt N_fc_iv]);

fc_iv0_dev2_f2 = reshape(permute(fc_iv0(:, dev2_ind(:, 2), :), [2 1 3]),...
    [N_dev2*N_mkt N_fc_iv]);

corr_g_ar = [gen_iv(dev_profit_all_dev2_f1, dev_profit_one_dev2_f1, ...
    profit_cutoff_corr, fc_iv0_dev2_f1), ...
    ...
    gen_iv(dev_profit_all_dev2_f2, dev_profit_one_dev2_f2, ...
    profit_cutoff_corr, fc_iv0_dev2_f2)]; 


corr_g_ar_boot = zeros(N_mkt*N_dev2, size(corr_g_ar, 2), N_boot);

for nb = 1 : N_boot

    
    corr_g_ar_boot(:, :, nb) = [gen_iv(dev_profit_all_dev2_f1_boot(:, nb),...
        dev_profit_one_dev2_f1_boot(:, nb), ...
        profit_cutoff_corr, fc_iv0_dev2_f1), ...
        ...
        gen_iv(dev_profit_all_dev2_f2_boot(:, nb),  ...
        dev_profit_one_dev2_f2_boot(:, nb), ...
        profit_cutoff_corr, fc_iv0_dev2_f2)];

end

%% get indices for the chosen sample
all_ind = reshape(1:N_mkt*N_dev, [N_dev N_mkt]);  
dev2_ind_f1_tr = all_ind(dev2_ind(:, 1), :); dev2_ind_f1_tr = dev2_ind_f1_tr(:);
dev2_ind_f2_tr = all_ind(dev2_ind(:, 2), :); dev2_ind_f2_tr = dev2_ind_f2_tr(:);

%% get the actual outcome of the pairs
entry_dev2 = entry_sim(:, dev2_ind(:, 1)) .* entry_sim(:, dev2_ind(:, 2));
entry_dev2 = entry_dev2'; entry_dev2 = entry_dev2(:);

%% market-level fc shock draws
corr_sim = reshape(...
    repmat(randn(1, N_mkt, N_corr_sim), [N_dev 1 1]), ...
    [N_dev*N_mkt, 1, 1, 1, N_corr_sim]);



